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Abstract 

A new Jacobian approximation is developed for use in quasi-Newton methods 
for solving systems of nonlinear equations. The new hypersecant Jacobian ap- 
proximation is intended for the special case where the evaluation of the functions 
whose roots are sought dominates the computation time, and additionally the 
Jacobian is sparse. One example of such a case is the solution of the discretized 
transport equation to calculate particle and energy fluxes in a fusion plasma. 
The hypersecant approximation of the Jacobian is calculated using function val- 
ues from previous Newton iterations, similarly to the Broyden approximation. 
Unlike Broyden, the hypersecant Jacobian converges to the finite-difference ap- 
proximation of the Jacobian. The calculation of the hypersecant Jacobian ele- 
ments requires solving small, dense linear systems, where the coefficient matrices 
can be ill-conditioned or even exactly singular. Singular-value decomposition 
(SVD) is therefore used. Convergence comparisons of the hypersecant method, 
the standard Broyden method, and the colored finite differencing of the PETSc 
SNES solver are presented. 

Key words: quasi-Newton method, Broyden's method, singular-value 
decomposition 



1. Introduction 

The Newton- Raphson method can provide solutions to systems of nonlinear 
equations. It works by finding a linear approximation to the residual vector (the 
vector whose root is desired), and then solving the resulting linear system to 
obtain an approximate solution. Iteration of this process, if it converges, gives 
a numerical solution to the nonlinear equations. Thus, there are three parts to 
each iteration. The first is to obtain the linear approximation, i.e., the Jacobian, 
the matrix of the derivatives of the residual vector components with respect to 
the parameters. The second is solving the resulting linear system. The third is 
to correct a change that makes the solution worse (e.g., by using a line search 
along the direction of the change). Here we concentrate on the first two parts. 
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Generally, calculating the Jacobian is the most computationally intensive 
part of the problem because the the Jacobian has TV 2 terms for N variables 
(and N residuals) . In many cases of interest it cannot be calculated directly, as 
the residuals are not known as simple analytic functions. Hence, the Jacobian 
is found by the finite difference approximation. However, when the number 
of variables is large, this becomes prohibitive, hence the development of quasi- 
Newton methods, such as the Broyden method [l[ . In the Broyden method, the 
Jacobian is updated by adding a term to it so that it accurately predicts the 
change in the residuals of the previous iteration. However, the Broyden method 
can be problematic; it is often found not to converge, in large part due to this 
update not converging to the true Jacobian. 

Our motivating application is the time advance of systems of nonlinear par- 
tial differential equations. In such cases there are a number of variables as- 
sociated with each cell upon discretization, and the residual is the change of 
those variables in a given time step that gives one a stable advance. Since the 
continuum equations are local, the residuals for a cell depend on the values of 
the variables in that cell and some nearby cells. Hence, the Jacobian is sparse 
with O(N) terms, and so it is reasonable to find the Jacobian through finite 
differences. In particular, one can use coloring, as is done in the PETSc SNES 
solver to ensure that the minimal number of evaluations are carried out to 
obtain those finite differences. 

We combine ideas from both these two approaches. We recognize that the 
Jacobian is sparse, and so there are only a small number of derivatives to be 
calculated for each residual. However, rather than explicitly computing all the 
finite differences at each step, we use the results from the last several solver iter- 
ations to solve for an approximation to the Jacobian. This may not be enough 
information for a precise solution, as either there may be an insufficient number 
of equations, or some of the equations may be degenerate (as can happen when 
successive steps are in the same direction). In the former case, which occurs for 
the first few iterations of a quasi-Newton solve, we use Broyden until sufficient 
state has been generated to calculate the hypersecant Jacobian. In the latter 
case we use singular-value decomposition (SVD) to determine the Jacobian. Fi- 
nally, because the resulting Jacobian may in general itself be ill-conditioned, 
we solve for the linear approximation for the change in the parameters again 
through SVD. 

We compare our method, the Broyden method, and the (colored) PETSc 
SNES solver to three problems. The first is a three- variable linear system, which 
illustrates how the Broyden Jacobian does not correctly converge. In contrast, 
our hypersecant method does provide the converged Jacobian. For linear prob- 
lems, once the Jacobian is correct, the problem is solved in the next step. Thus, 
for this case there is no need to apply the PETSc SNES solver as the results 
are known. The second problem is a three-variable nonlinear system. There 
we find that the Broyden and hypersecant methods are competitive, with the 
finite-difference approximation taking more function evaluations. This is per- 
haps not surprising, as with only three variables, the Jacobian is nearly dense. 
Finally we apply the methods to our motivating application, a nonlinear trans- 
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port equation of a kind that arises in certain problems in plasma physics. Here 
we find that the hypersecant method takes roughly half the number of function 
evaluations to converge compared with the finite difference approximation. The 
Broyden convergence is much slower still. 

The organization of this paper is as follows. In the next section we discuss 
the Broyden method for updating the Jacobian. In the subsequent section we 
develop our hypersecant method. Next we compare the methods in application 
to the simple three- variable linear and nonlinear problems. Subsequently we 
compare them in application to our motivating problem, a discretizcd transport 
equation. Finally we summarize. 

2. Quasi-Newton methods and the Broyden approximation 

The goal of any Newton-like method is to solve a system of N nonlinear 
equations 

F(x) = 

of N variables, x. Given a current best guess x fc for the solution, one Taylor 
expands in its vicinity to obtain the linear equation, 

F(x fc + <5x fe ) = F(x fc ) + J(x fe ) • <5x fe + 0(6x k • <5x fc ) = , (1) 

where J is the Jacobian 

dF 

which implies that a better approximation x k+1 to the root is found by solving 
for the linearized change 5x. k , 

J(x fc ) -Sx k = -F(x fc ), 

from which one obtains the improved guess for the root, 

x fe+i = x fe + Sx k 

Iterating this process then gives the root, if the process converges. 

In the standard Newton-Raphson method the true, analytic Jacobian is used. 
If the true Jacobian is not known, a finite-difference approximation can be used, 
making the scheme a quasi-Newton method. However, if F(x) is numerically 
expensive to evaluate, the finite-difference Jacobian can become inefficient due 
to the number of function evaluations needed to compute it. 

The Broyden method [l[ is a common quasi-Newton method where the Ja- 
cobian is approximated as: 

Jfc+1 = J" + A A ® . ( 2 ) 

where we have introduced the shorthand notation F fc = F(x fc ) and J fe = J(x' c ), 
and where £g> is the dyadic (or outer) product. The above approximation comes 
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from demanding that the new Jacobian J fe+1 predict the change in the residuals 
for the most recent step, 

F fc+i _ F fe = jfe+i . (- x fc+i _ x fc) _ 

The initial approximation J° is typically set to either the finite-difference Jaco- 
bian (good, but numerically expensive) or the unit matrix (free). 

The obvious advantage of the Broyden Jacobian is that it is free in the sense 
that it requires no extra evaluations of F. The disadvantage of the Broyden 
Jacobian is that it is not guaranteed to converge toward the true Jacobian; in 
fact error can accumulate until the increment Sjc found by inverting it is not 
even in the descent direction, i.e. x + <5x is not closer to the root than x. This 
occurs because Eq. @ is underdetermined, so the update ^ is only one of many 
possible. In a solver using the Broyden approximation one must intercept this 
occurrence and reinitialize the Jacobian. Indeed, in our experience from using 
the Broyden method to solve a discretized fusion transport equation Q , optimal 
performance, i.e. the minimal number of function evaluations, required the 
Broyden Jacobian to be reinitialized as often as every 5-10 Newton iterations. 

3. The hypersecant approximation of the Jacobian 

Our hypersecant approximation is like the Broyden method in that the Ja- 
cobian comes with no additional function evaluations, but unlike Broyden, the 
approximate Jacobian will approach the true Jacobian after sufficiently many 
iterations as one approaches the solution, provided the directions of the steps 
span the relevant space. The primary difference is that instead of requiring only 
[21 we require that this hold for several previous steps, 

F(x fe+1 ) - F(x fc ~0 = J fc+1 -(x fc+1 -x fe ^), t = Q,...,L-l. (4) 

For the fully dense Jacobian, one obtains a full set of equations after N linear 
independent steps. For large systems, this is not practical in general. However, 
it can be very good when the Jacobian is sparse, such as occurs in the implicit 
advance of discretized partial differential equations. In particular, the above 
system of equations is determining when the sparsity is such that any residual 
depends on only L variables. For example, in a finite difference discretization of 
a one-dimensional, single- variable, second-order PDE, each residual depends at 
most on three variables, and so after three linearly independent steps one has a 
good approximation to the Jacobian. 

One issue is that the above approach cannot alone be used to find a new 
Jacobian via Eq. ((4]) until one has taken a number of steps equal to the maximum 
number of nonzero elements in a row of the Jacobian. This is clearly the case 
for the first few iterations. The Broyden update can then be used to reduce the 
number of unknowns to make the system well determined, see Section |4~T1 until 
a sufficient number of steps have been taken. Even if a sufficient number of 
steps have been taken, the system of equations ((4]) can be underdetermined, 
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e.g. because several successive steps are in nearly the same direction. This 
case is handled by using singular-value decomposition (SVD) to solve for the 
hypersecant Jacobian elements, see Section ETT1 



4. Comparison of hypersecant, Broyden and PETSc colored FD Ja- 
cobian for few-variable systems 

In this section we compare the application of two approaches (hypersecant 
and Broyden) to a linear three-variable system, and then we compare the ap- 
plication of three approaches (hypersecant, Broyden, and finite-difference Ja- 
cobian) to a nonlinear three-variable system. Our solvers did not attempt to 
correct for an increment that takes the solution further away from the true root. 



4-1- Linear system 

The linear system is defined by the equations, 



h(xi,x 2 ) 
f 2 {x 1 ,x 2 ,x 3 ) 
/$( x 2 , x 3 ) 
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This system has the true solution x = (xi, x 2 ,x 3 ) 
is 

1 1/2 



1/2 




(1,1,1), and the Jacobian 



(5) 



1 1/2 
1/2 1 

A simple quasi-Newton solver was implemented using the LAPACK Q im- 
plementation ATLAS Q version 3.8.3. This simple linear system could of course 
have been solved in a more direct manner, but the purpose of the test is pri- 
marily to compare the Broyden and hypersecant Jacobian approximations. 

We initialize the Jacobian as the unit matrix, J = 1 and use the initial 
guess x = (1/2,1/2,1/2). It takes 14 Newton iterations using the Broyden 
approximate Jacobian to calculate the root to 12 significant digits. The reason 
for the poor convergence is that the Broyden Jacobian does not converge toward 
the true Jacobian. At the final iteration, the Broyden Jacobian is 



B 



fc=14 



1.124366 0.599240 
0.397825 0.823632 0.397825 
0.599240 1.124366 



which is far from the exact Jacobian (O . 

For the hypersecant solve, we initialized the Jacobian with the identity and 
used the Broyden update for the first two iterations, but after the third iteration 
we used Eqs. ((4]). As an example, for the middle row, the system of equations 
to solve is: 
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Thus, at the third step we obtain the hypersecant Jacobian, 



H 



fc=3 



1.000000 0.500000 
0.500000 1.000000 0.500000 
0.500000 1.000000 



It is seen to equal the finite-difference approximation of the Jacobian (which 
in turn equals the true Jacobian for the linear system used in this case). This 
simple example thus clearly illustrates how the hypersecant Jacobian avoids 
the major disadvantage of the Broyden approximation, while maintaining the 
greatest advantage: the fact that it is calculated without making any extra 
evaluations of the function f (x) . 

This linear example also illustrates how singular-value decomposition (SVD) 
is used to solve for the hypersecant Jacobian elements. The variables x\ and 
x-i are interchangeable. So if they are given the same initial values, they will 
follow exactly the same trajectory. The linear system to solve, Eqs. (J4j) for the 
hypersecant Jacobian elements on row 1, therefore becomes singular: 



-1.429236794928 x 10" 
-3.617952748235 x 10" 
+6.382047251765 x 10" 



-4.143137616670 x 10" 
-3.254780687737 x 10" 
+4.245219312263 x 10" 



Hi 
Hi 



-4.143137616670 x 10" 
-3.254780687737 x 10" 
+4.245219312263 x 10" 

-1.843550556595 x 10" 1 
-6.872733435972 x 10 _1 
+ 1.062726656403 x 10 +0 



Column 1 and 3 are identical and LAPACK dgesvxO fails with error code INFO 
= 4 signifying a coefficient matrix singular to within machine precision. The 
inverse condition number RC0ND = 7 . 426728596992e-17 confirms this. The 
implementation of the hypersecant solver used here handles situations like this 
by switching to SVD when dgesvxO fails. LAPACK dgesvdO informs us that 
the singular values of the coefficient matrix are 



(wx,w 2 ,w 3 ) = (1.060808064514, 9.516998001847xl0" 2 , 4.162068222172x10 



-17\ 



Following the usual SVD prescription to solve linear systems, our hypersecant 
solver implementation treats W3 as approximately zero and sets 1 /1V3 = when 
it calculates the hypersecant Jacobian elements 



/ Tj3 Tj3 Tj3 

1^1.1^1.2^1. 



(0.500000000000, 1.000000000000, 0.5000000000000) 



which equals the true Jacobian to within machine precision, exactly the answer 
we want. The fact that Eqs. can be ill-conditioned or even exactly singular 
therefore does not mean that their solutions, the hypersecant Jacobian elements, 
are a poor approximation. As shown here, by using SVD, hypersecant Jacobian 
elements can be extracted, exact to within machine precision, even from a very 
ill-conditioned system. 
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4-2. Nonlinear system 

The next test is for the system of nonlinear equations 



k 

4 



f2(xi,X 2 ,X 3 ) 
/»( X 2 ,Xz) = 

This system also has the root x 
the root is 



(7) 



(x±, X2, Xa) = (1,1,1), and the Jacobian at 



1 1/2 
J(l,l,l)= 1/2 1 1/2 
1/2 1 

To be able to compare not only the hypersecant and Broyden Jacobians, but 
also a finite-difference (FD) Jacobian, a quasi-Newton solver was implemented 
using the PETSc SNES nonlinear solver toolkit. PETSc can do colored finite 
differencing, i.e. take full advantage of the sparsity of the Jacobian to minimize 
the number of function evaluations used to calculate the FD Jacobian. We only 
specify the sparsity pattern of the Jacobian and let PETSc determine how to 
do the colored finite differencing. 

For the hypersecant Jacobian, we take advantage of the fact that the first 
and last rows of the Jacobian have only two non-zero elements. For the first 
row, the hypersecant Jacobian elements are given as the solution to the linear 
system 
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For the first and the last rows of the Jacobian, the hypersecant approximation 
can thus be used already after two iterations (the hypersecant Jacobian elements 
can be calculated from the available state at the end of the second iteration, 
and is thus first available for use during the third iteration). 

When the iteration number is too small to calculate the hypersecant Jacobian 
elements (k = 1 for row 1 and 3, and k — 1, 2 for row 2), we first do a Broyden 
update of the old Jacobian. Elements of this updated Jacobian are then used to 
replace sufficiently many unknowns in Eq. ((4J to reduce the number of unknowns 
to equal the number of iterations, thereby making the system well determined. 
For example, for the first iteration we get for the first row of the Jacobian 

HlM- x l)+ B lM-xi) = fi-f l , 

which only has the single unknown H\ x , where the tilde is used to indicate that 
this is only an approximate hypersecant Jacobian element. This equation can 
be trivially solved, 
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Figure 1: Accuracy as a function of the number of function evaluations when doing quasi- 
Newton solve of Eq. J7J. For the triangles, the Broyden approximation of the Jacobian was 
used. For the diamonds, the hypersecant Jacobian. For the squares, the PETSc colored 
finite-difference Jacobian. 



Because the diagonal elements have larger values, we somewhat arbitrarily 
choose to use B\ 2 to calculate H\ x , instead of B\ A to calculate H\ 2 - Sim- 
ilarly, for the middle row we get 

tt! _ ~ fi ~ ^2,1 ( x l — x l) ~ ^2,3 ( x 3 ~ x 3) 

etc. For k > 3, the full hypersecant Jacobian is calculated (and subsequently 
used for k > 4). 

As usual, wc initialize the Jacobian as the unit matrix and use the initial 
guess x = (1/2,1/2,3/2). As can be seen in Fig. [1] the quasi-Newton solve 
converges to the default SNES accuracy in eleven iterations when either the 
hypersecant or Broyden Jacobians are used. Using the SNES colored FD ap- 
proximation of the Jacobian, the solve converges in just five iterations. However, 
because the colored FD approximation requires three extra function evaluations 
(for a Jacobian with this particular sparsity pattern), the total number of func- 
tion evaluations is 21 vs. the 11 needed by hypersecant/Broyden. 

Even if the final Broyden Jacobian is a much poorer approximation of the 
true Jacobian than the final hypersecant Jacobian is, the performance of the 
Broyden and hypersecant methods in this case is very similar. However, the 
hypersecant Jacobian gives more rapid convergence for the last few iterations, 
as it becomes a more accurate approximation to the true Jacobian. 



8 



5. Application to a nonlinear transport equation 



As noted in the introduction, our motivation for this research is to develop 
methods for solving nonlinear transport equations. Our particular example 
comes from cross-flux-surface transport in a fusion confinement device, where 
the transport equation is a radial continuity equation [f| : 

where the highly nonlinear flux is given by T = T(r, u, u') and u' = du/dr. The 
volume clement, the differential volume in the radial variable r, is denoted by 
V'(r). The components of u(r) are typically the temperatures and densities of 
the plasma particle species, but they could also include momenta. In the cases 
where it is difficult to advance this partial differential equation, the dominant 
contribution to the flux comes from turbulent transport, which is triggered 
when a critical gradient is exceeded. This creates a sensitive dependence on u'. 
The stiffness of the transport equation makes it necessary to use implicit time 
discretization. 

Time-implicit, spatially second-order finite-difference discretization on an 
equidistant mesh Tj = j Ar, j = 0, . . . , N gives the discretized system of nonlin- 
ear equations: 

Fj-(Auj_i, Auj, Au i+ i) = Auj + 

J i ^^ + i/ 2 ) r ^/ 2 -^(^-i/ 2 ) r ;-i/ 2 cn+1 x 



- S" +1 At- 



\W'(rj) Ar > J (io) 

( 1 W , (r j+1/2 )T r },-,^-W'(r i _ 1/2 )T^ 1/a \ 

(!-*) T J ; t - ' J ' 1/2 -S" At = 0, 

\W'(rj) Ar 3 j 

where subscript denotes spatial location and superscript time level, e.g. S™ = 
S(fj,nAt) is the source at radius r = Tj and time t = nAt. The indepen- 
dent variables Au 3 = u™ +1 — u" and T J+1 / 2 = r(r J+1/ / 2 , u J+1 / 2 , u^. +1 ^ 2 ), with 
u i+i/2 = ( u j + u i+i)/2 and u'- +1 ^ = (u J+ i — Uj)/Ar. The implicitness pa- 
rameter 6 — gives fully explicit (FE) equations, 9=1 fully implicit (FI) and 
6 = 1/2 time-centered implicit (CI), with second-order temporal accuracy. 

We next compare the three Jacobian approximations for this transport equa- 
tion (fT0|) with a single profile evolved [e.g. the ion temperature 7i(r)]. For the 
nonlinear flux we use the linear critical gradient model, which assumes the flux 
r = — \u' and diffusivity 

X = max{(l/L - 1/L C )/L, Xmm} , 

where L = u/\u'\ and L c and Xmin are parameters. We also use a source 
term S(r) = 1 — r a . We somewhat arbitrarily choose the parameters L = 1/2, 
Xmin = 1/10 and a = 2. With fully implicit time discretization we get the 
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system of nonlinear equations 



Fj(Av,j-i, Av,j, AMj+i) = Auj+ 

( 1 W>(r J+1/2 )r^ /2 - W'jr^r^ \ (11) 

[w^~) At % ) At ~°> 

j = 1, ...,N — 1. To determine the appropriate on-axis (ro = 0) boundary 
condition, we take one step back and write the flux term as an undiscretized 
divergence term: 

W (r )F = W'(r )Au + (JLw'(r)T) At - W\r Q )S^ +1 At = , 

where we have also multiplied the equation with W'iro). In the limit ro — > 0, 
in which W'(ro) — > 0, we get 

d \ 
T-W'(r)r At = 0. 
dr J Q 

Discretizing this equation with second-order error gives 

(12) 

which is a nonlinear on-axis boundary condition. The first term depends on Aito 
and Aui, and the second term on Ami and Au 2 - The first row of the Jacobian 
will thus have non-zero elements in the first three columns. The Jacobian is 
then tri-diagonal, but with Jo, 2 ^ 0, due to the nonlinear BC. 

The hypersecant Jacobian is calculated in the way described above in Sec- 
tion [4721 but both hypersecant and Broyden Jacobians are initialized as the unit 
matrix with the first row replaced by Jo.o = 3/10 x N, Jo,i = —4/10 x N and 
J0.2 = 1/10 X N (the values given by the nonlinear BC at steady state). 

Taking a single time step of 0.1 ms, the PETSc SNES quasi-Newton solve 
converges as shown in Fig. [5] for the three different Jacobian approximations. 
With the hypersecant Jacobian, the solve converges using 9 function evaluations. 
With PETSc colored FD 16, and with Broyden 20 function evaluations. With 
PETSc colored FD, only three quasi-Newton iterations are needed, but during 
each iteration, four extra functions calls are used to calculate the colored-FD 
Jacobian. With hypersecant Jacobian, the convergence is poor until the fourth 
iteration, the first iteration where the full hypersecant Jacobian is used. After 
that it takes five more iterations until convergence (for a total of nine), but each 
iteration only requires a single function evaluation. With the Broyden Jacobian 
the convergence is much poorer. 

For larger time steps, a problem with using the hypersecant Jacobian ap- 
proximation for PETSc SNES quasi-Newton solves becomes evident. The hy- 
persecant Jacobian is a poor approximation of the true Jacobian for the first few 
iterations. To use hypersecant quasi-Newton solves for the transport equation 
with large time steps, one probably even needs to accept that the accuracy gets 
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Figure 2: Accuracy as a function of the number of function evaluations when doing quasi- 
Newton solve of Eq. lilt . For the triangles, the Broyden approximation of the Jacobian 
was used. For the diamonds, the hypersecant Jacobian. For the squares, the PETSc colored 
finite-difference Jacobian. 



worse for the first few iteration. The purpose of the first few iterations is not 
primarily to get closer to the root, but to sample enough function values to 
make Eq. Q well determined to allow the calculation of the full hypersecant 
Jacobian. PETSc SNES tries to improve the accuracy for every iteration by do- 
ing a line search for the minimum accuracy along the direction of the increment. 
This typically fails for a poor approximation of the Jacobian. We are currently 
investigating some promising ways to solve these problems and plan to report 
on our findings in a future publication. 



6. Conclusions 

We have introduced the hypersecant Jacobian approximation, which has im- 
proved characteristics for solving systems of nonlinear equations with sparse 
Jacobians, such as occurs in the implicit advance of nonlinear PDEs. The basic 
idea behind this approximation is to keep several historical values - sufficient to 
be able to solve uniquely for an approximate Jacobian generically, but with use 
of Singular Value Decomposition to find a best solution when the equations are 
degenerate. Like the Broyden Jacobian approximation, the hypersecant Jaco- 
bian is calculated without making any additional function evaluations. Unlike 
the Broyden Jacobian it converges toward the finite-difference approximation of 
the Jacobian. 
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Multiple examples (linear and nonlinear systems of three independent vari- 
ables and the ID transport equation with linear-critical-gradicnt diffusivity) 
were solved using these different methods. The new hypersecant method was 
shown to lead to a convergent Jacobian. Improved overall convergence with 
respect to number of function evaluations was observed. 
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